# Dickstein, Ho, and Mark (2023)
# This file defines the cost variables needed in the premium setting equation for each bootstrap draw

# Preliminaries
setwd("../../../library")
source("PreliminariesCode.R")
library(readr)
library(knitr)
library(readxl)
library(numDeriv)
library(yaml)
library(haven)
library(Hmisc)
library(EnvStats)

ravec <- 1:7
yrvec <- 2014:2016
k_vec <- 1:20

RNGkind("L'Ecuyer-CMRG")
generate_randomseed <- function(run){
  s <<- .Random.seed
  if(run > 1){
    for (i in 1:(run - 1)) {
      s <<- nextRNGStream(s)
    }
  }
  .Random.seed <<- s
}

before_stuff <- c(ls(), "before_stuff")

# Start loop
for(k in k_vec){
  
  # In loop preliminaries
  solution_path <- paste0(project_folder, "/analysis/estimationcode/specs/lowrisk_simplemh_censor0.2_main/bootstrap/output_", k) # Where the demand system solution lives
  output_path <- paste0(project_folder, "/analysis/prem_setting/bootstrap/output") # where the results are saved
  
  # Load the necesary data
  
  ## Demand system solution
  solution_csv <- read_csv(paste0(solution_path,  "/solution.csv"))
  covar_names_csv <- read_csv(paste0(solution_path, "/covar_names.csv"))
  covar_lens_csv <- read_csv(paste0(solution_path,  "/covar_lens.csv"))
  solution <- as.numeric(solution_csv[,3:ncol(solution_csv)][1,])
  covar_names <- covar_names_csv[['x']]
  covar_lens <- covar_lens_csv[['x']]
  
  ## Exploded data
  source(paste0(project_folder, "/library/DA01QuanModifyExplodedDFunc.R"))
  explodeddata <- fread(paste0(project_folder, "/data/explodeddata.csv"))
  
  ### Reconstruct to the exploded data we used for this bootstrap draw
  set.seed(1)
  generate_randomseed(k)
  samp <- sample(unique(explodeddata$hhid), length(unique(explodeddata$hhid)), replace = TRUE)
  setkey(explodeddata, "hhid")
  explodeddata_k <- explodeddata[J(samp), allow.cartesian = TRUE]
  ### Create modified exploded data
  acgquinsvec <- as.numeric(scan(file=paste0(project_folder, "/analysis/orig/acg_positive_quintiles"), sep = ",", what = "character")[2:7])
  print(paste0("Modifying... ravec: ", paste(ravec, collapse = " and "), "yrvec: ", paste(yrvec, collapse = " and ")))
  data <- mod_exploded_data(
    explodeddata_k, 
    paste0(project_folder, "/data"),
    highmeanthresh = 0.59616,
    lowsumthresh = 0.2575,
    highsumthresh = 1.76329,
    highmaxthresh = 3.37157,
    medincome = 2.5,
    acgquins = acgquinsvec,
    winsorbound = 89.3707340916662,
    ravec = ravec,
    yrvec = yrvec)
  data <- data.frame(data)
  
  ### turn this into the "parametrized cost data" format
  source(paste0(project_folder, "/analysis/estimationcode/parametrize_choice_data.R"))
  param_data <- parametrize_choice_data(solution,
    data,
    covar_names,
    covar_lens,
    score=TRUE)
    
  # Draw from lambda distribution and simulate costs for plans [from generatestderrors.R]
  set.seed(4197169)
  param_data <- param_data %>%
    mutate(lambda = rexp(n(), rate=alpha_i),
      simcost=(x_av+omega_i*x_av^2)*lambda,
      wsimcost=simcost*(simcost < 89.37073) + 89.37073*(simcost > 89.37073),
      wlambda=lambda*(lambda < 89.37073) + 89.37073*(lambda > 89.37073))
  param_data['sim_eps'] <- revd(nrow(param_data))
  param_data['sim_util'] <- param_data['delta_i'] + param_data['sim_eps']
  param_data <- param_data %>%
    group_by(subscriberid_year) %>%
    mutate(sim_choice = as.numeric(sim_util == max(sim_util))) %>%
    ungroup()
  
  # Collapse and create plan cost data
  source(paste0(project_folder, "/analysis/prem_setting/bootstrap/predict_cost_fun.R"))
  predict_cost(seed = 4197169, 
    param_data = param_data, 
    output_path = output_path, 
    k0 = k)
    
  # Cleaning up before next run
  rm(list=setdiff(ls(), before_stuff))
  print(paste0("This is in the environment:", ls()))
# End Loop
}
